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SYNOPSIS 


A knowledge of the way in which a flood wave 
attenuates as it travels downstream along a long open 
channel with a continuously moving unsteady- steady 
boundary is of importance in flood routing. 

Implicit scheme was used and a dimensional 
analysis was carried out to consider all the significant 
parameters that affect the flood propagation, A systematic 
analysis was made by varying all the parameters in the 
practically possible range. 

It was found that relative flood wave amplitude 
depends upon initial wave amplitude ratio of time to peak 
to time to centre of gravity of the inflow hydrograph, 
and Manning's n. The first two terms take into account 
the rate of rise of hydrograph, time of travel of the 
wave peak is chiefly dependent upon initial wave amplitude. 
These conclusions are very important in prediction of the 
peak stage and time of its occurrence. 
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CHAPTER I 
ITOODTJCTION 

Our civilization has grown on the hanks of rivers 
and for hundreds of years they were and a,re still a 
boon to mankind. But with the development of science, 
man tried to bind these rivers to flow on a particular 
path, wanted to cross them through strong bridges and 
thus to a great extent disturbed their natural flow. 
Some times dams are constructed, water is regulated in 
a particular way and thus taking the control in one's 
hand. One of the most important of all the unsteady 
flow phenomena the engineer has to deal with, is the 
movement of a flood wave down a channel, usually in 
the case of a river, and the problem along with this 
is the tracing of this movement and any related changes 
in the form and height of the wave. It is necessary 
that the engineer should possess theoretical or 
accurate semi-empirical means of determining the 
behaviour of a flood wave in a channel of specified 
form and slope, so that he may be able to predict 
the effect of the changes v^hich he ofcen makes in a 
channel in the interests of channel improvement and 


flood control r 
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Flood routing is the process of tracing hy 
calculation the course of a flood wave and in large 
measure this problem is an application of the principles 
of unsteady floy;. However a flood wave has its own 
features and may req.uire its own special techniques for 
exploring these features. As for example, rise and 
fall of a flood occurs much more slov\rly than many of 
the flow changes in unsteady flow phenomena and so some 
of the acceleration terms in the equation of motion may 
not be of much significance, or otherwise, when the 
flood wavs propagates, it becomes longer and lower 
as it moves downstream and this subsidence is very impor- 
tant for a practicising engineer. He may, therefore 
prefer instead of a general solution, an approach 
which concentrates on those aspects of the solution 
that are particularly relevant to the subsidence problem. 

A semi-empirical method known as the Muskingum 
method, first developed by McCarthy [3I ] has been in 
use for many years. This method, which uses the 
equation of continuity and some empirical constants 
which characterise the basin concerned y/as presented by 
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Linsley [20| . This coefficient method is by no 
means an exact method of routing hydrographs in as 
much as only a few of terms of governing equations are 
considered in the derivation. The method though adequate 
for many purposes such as in the planning and designing 
stages of flood control or multiplex-purpose projects, 
reduction in natural flows resulting from the control 
of floods by reservoirs etc., is of limited utility as 
it d0'-'’s not take into account the dynamics of flow. 

In recent years the advancement of computer technology 
has estimulated the development and application of 
a more rigorous approach of the numerical integration 
of the equations of motion and continuity for unsteady 
flow. A great deal of effort has been made to develop 
computational methods applicable to the routing of a 
flood wa.ve in a complex natural or hypothetical drainage 
system. However, it is worth mentioning that the 
significant factors influencing the unsteady flow 
phenomena are so numerous and any typical watershed 
is so complicated that it is very difficult to 
arrive at any reliahle conclusion. 
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The precise role of each factor a.ffecting the 
flood wave propagation in open channels can best be 
understood by studying the flow in a relatively simple 
and ideal channel in ¥/hich all important parameters 
can be independently and systematically vaxied. 

For example, the relations a,mong important flood 
cha,racteristics, such as maximum stage vjith Ifenning’ s 
roughness, initial wave amplitudes .time to peak of the 
hydrogreiph, distance etc. may be studied. It may be 
also possible to a.rrive at a result which can be used 
to get a„n empirical formula and thus help in flood 
prediction. 

An approach of this type in open channel was 
applied by I/P.zayaney and Song. But their work is of 
very restrictive nature as the wave applitude and, the 
length of the channel are small. 

The transient flows due to a +v(ly skewed 
hydrograph of given duration at the upstream end of 
long rectangular straight channel were computed by 
implicit method. Implicit method was used because it 
is a very stable method and can be used very easily in 
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the case of nature,! rivers as shown hy Amein In his 
paper | 4 1 • Though the stability of the scheme makes 
its use convenient and easy for various dT and dx i.e. 
time step and distance step, it will have some 
numer i c al damp ing . 

The inflov; hydrograph taken was +vely skewed 
one because most of the practical hydrographs resemble 
this shape. Sine curve can also be used but then 
it does not represent the actual field problem and 
since the motive is to study subsidence with more 
practiccil input data, log poarson type III hydro graph 
was considered more useful. The initial wave amplitudes 
were 0.5, 1 and 2 times the initial water depth, 
the ratio of time to peak to the time to center of 
gravity of the hydro graph was l/3 and 2/3, 

Manning's roughness coefficients were 0.02, 0.025, 0,03. 
These parameters were used in different groups and 
the resulting propagation of the flood v/ave was ' 
calculated. The damping of the peak stage and the 
time of o.ccurrence of peak v/ere studied in detail 
which will be presented in this report. 



CHAPTER II 



2, 1 Assumptions 

The one dimensional eq_uations of continuity 
and momentum of unsteady flow in open channels are 
referred to as St. Tenant ’s equations. These equations 
are derived in any standard text book on unsteady flow 
e.g. Mahmood and Yevjevich, Chow, Henderson etc. The 
following are important assumptions in the derivation 
of these equations, 

( 1 ) The channel is straight and uniform in the 
reach so that the one dimensional approach can 
be applied. 

(2) Telocity distribution across the wetted area 

does not substantially affect the wave propagation. 

(3) Priction losses are given by Manning’s formula 
i.e. it is equal to loss in steady flow conditions, 

(4) The average slope of the channel bottom is 
small therefore tan a = Sin a and Cos a = 1 , 

(5) The wave surface varies gradually which is 
similar to stating that the pressure distribution 
is hydrostatic. 
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(6) E1oy>; is suTDcritical so the kinematic wave is 
important. 

Sq.uation of Oontinuity; 

A short reach of channel of length x with 
the flow taking place from section 1-1 to section 2-2, 
shown in figure 2,2,1, Assuming x to represent the 
horizontal distance, p the density of water, g, the 
acceleration due to gravity, z the channel bottom 
elevation, y, the water depth, A the cross-sectional area; 
and V the average velocity. If Q be the volume flow 
rate entering the channel through section 1-1; 

Q + ( f^x will be the volume outflow rate; the 

mass f water entering the channel reach during 
time will be p- Q. £it. 

The mass of water leaving the channel during 
the same time would be 

P (Q + ^ x) ■ h-t. 

Assuming Q to increase with xi,e. 6Q/8x is +ve 
there is a net mass outflow from the reach. Applying 
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law of mass conservation during time At, it can be 
easily concluded that mass inside the reach will be 
reduced and consequently the water surface should fall. 
The mass of water inside the reach of length Ax is 

X A Ax. 

The rate of decrease of this mass is 

- ~ (fAAx) = - Pax~ 

dt 9t 

The reduction in mass during time interval t is 

^ P C, X ( ?A-) £X t 

0t 

for the conservation of mass, 

the net mass outflow = Reduction in mass inside the 

P ( Q + h x) At'-('? QAt) = .-p At AX 
which gives on simplification 

^ 
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Since Q = AV 


Y 


6A 

5'x 


+ 


A 


ay . 8A 
ax at 


0 


(1) 


The inomentura equation can "be easily derived from the 
Newton’s second law of motion which states that the 
time rate of change of momentum is equal to the external 
force applied, both being considered in the some direction. 


Referring again to Fig. 2,1, the applied forces 
on the element are pressure force, gravity and frictional 
forces. If water depth is y at section 1-1, water 

0 V 

depth, at 2»-2 v^ill Y ^ ^ ^ simila.rly area 

at 1 -1 is A, area at 2-2 is A + ^ Ax assuming 
higher order terms f Taylor’s series vanish or are 
negligibly small. 

Pressure at section 1 - 1 acting to right = P g . y 


and the same at 2-2 acting to left, 

= f g ( y + “ Ax) 

8x 


Pgj - pg(y I" ^ 


‘g -21 X. 


net pressure force 


■f 
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Unbalanced pressure force is 

= - P g.i Di. X 

6x 

This will act towards the negative x-axis. 

The component of the gravity force which is in the 
direction of motion is 


= P gA hi X 


.6 z 
X 


Z 

but ■;-; — . = channel bottom slooo So, 

hx 

Therefore, the gravity force in the direction of motion 
is p gA Sq ih X. The shear force opposes the motion 
and is directed tovi^ard the left. 


Shear force can be expressed in terms of the head 
loss as follows, 

ZT P Ax = ^ ^ 

Where P is wetted perimeter and h^ is the head 
loss over distance A x. 

lij- = Sf Ax. 


So , ZT- P Ax 


Y AxA. 
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The resultant force on the element of volume in the 
direction of motion is the resultant of the graistational, 
frictional, and pressure forces: 

P SASq >- giSf X - f gA I"! .N, X 

The tot 8.1 rate of change of momentum is the sum of the 
local and convective momentum changes. 

The acceleration of the fluid element is sum of 
local and convertive acceleration 

§1 = n. + 

dt 6t 8x 

Therefore the product of mass 9 A and acceleration 

will give the rate of change -of momentum i.e. 


P A £ix 


{11 + Y ) = f g,is X 

et 8x 

- p gAS^ Cxx ^ 9 gil 5. X 

5x 


9Y 

ST 


+ 



So - 


h 

0X 


g S. 


( 2 ) 
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Equations (1) and (2) are the equations of unsteady 
flow. They serve as the basis of constructing the 
mathematical models for the investigation of numerous 
problems in unsteady flow. 

^ 5 Solu ti on of Governing Equat ions : 

2.3.1 Oompl ete Solutions; 

There are 5 approaches for solving the governing 
equations 

differential ^1) analytical solution (2) approximate 
numerical solution (3) complete numerical solution. 

It is only recently v/ith the advent of high speed 
computers it became possible to get the complete numerical 
solution in a comparatively less time and with greater 
accuracy, 

Numerical methods for the solution of differentia,! 
equations can be established by replacing the derivatives 
by finite differences. The differential equations 
are then represented by corresponding algebraic 
equations. In spite of the fact that numerical methods 
donot require any drastic simplifications of the 
equations of unsteady flow, the procedure for represent*^ 
ing the partial differential equations by algebraic 
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equations is not straight forward. Many factors 
involving convergence, stability, accuracy and efficiency 
come into picture. Although the procedure involved may 
be reduced to simple mathematical operations, the number 
of those operations may be extremely large, and this 
is true for application to most flood routing problems. 

Methods for the numerical, solution of the 
complete equations of unsteady flov-r based on finite 
difference schemes are broadly classified as direct 
methods and characteristic methods. Direct methods are 
those in which the finite difference representation is 
based directly on' the primary equations. In the 
characteristic methods, the equations are first trans- 
formed in to the characteristic from and the latter form 
is used in finite difference representation. A fixed 
mesh of points on the time-distance plane is employed 
to identify the points at which solutions are obtained 
by the direct methods. In the characteristic methods, 
solutions may be obtained at the intersection of 
characteristic curves on the time -distance plane or 
at fixed points of rectangular mesh by interpolation 
between the intersection of characteristics. 
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Tlie finite difference schemes used in the direct 
and characteristic methods are further divided into 
two - i.e. explicit and implicit methods. In the 
explicit methods, the finite -difference equations from 
which the unknowns can he evaluated explicitly a few 
at a time. The unknowns occur implicitly in the finite 
difference equations of implicit methods, v/hich are 
usually non-linear algebraic equations. 


Methods of Solution 


Direct 


Characteristic 


- j 

Explicit Implicit 




Characteristic 

nodes 


Explicit 


Implicit 


Fixed mesh 

I 

1 


Explicit Implicit 


Out of all these methods, only three are more commonly 
used. 

1. the implicit characteristic method using 
characteristic network. 
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2, the explicit direct method. 

3. the implicit direct method. 

The other methods have been found to be either unstable 
for general flood routing problems or they were basically 
sirailar to one of the three methods discussed above. 

2,3,2 I mplicit Method 

In the implicit method, x derivatives are 
replaced in terms of finite differences evaluated at 
time to + and thus the unknowns occur implicitly in 
the resulting algebraic equations, solution of these 
algebraic equa,tion is more complex as compared to the 
explicit method. 

There are many schemes which have been 
reported in the literature. Of these the following 
have been used for studying the open channel transients; 

iriessman’s scheme [_25j , .Amein’s scheme [4J , 

Vasilev's scheme [.26j etc. 

While writing the expression for partial 
derivative, a weighting factor a is introduced in the 
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Priessmann’s and in Amein’s schemes. The presence of 
a introduces artificial damping in addition to the 
damping due to friction and other losses. To make the 
scheme stable, value of a must lie between 0.5 to 1. 
Vasiliev's scheme has two steps and thus is more compli- 
cated and requires more computer time. Keeping in 
view the simplicity of Amein’s scheme, it was preferred 
to other ones. 

MP.zayaney and Song have used characteristic 
method. But in this study implicit scheme was preferred. 
The reasons are following: 

(1) There is no restriction on iii. t in this method. 
Since hydrographs with different rate of rise 

were attempted it was necessary to have some method 
which does not depend much on the nature of input 
data. 

(2) Economy: Because the size of fit is not 
stricted by any stability criterion, computer 
space requirement is less. 
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2,3.3 P etails of implicit Scheme 

On the (x,t) plane, the net points are determined 
hy the intersection of straight lines pa.rallel to the 
X and t axes. The lines drawn parallel to the t >- axis 
represent locations along the channel while those 
drawn parallel to the x-axis represent time. The 
spacing of locations has "been kept constant for the 
present analysis. The t axis is the upstream channel 
"boundary. The downstream "boundary condition is steady 
flow. 

The net work is drawn in Pig. 2,1, With reference 
to the Pig. 2,1, assumed that all the varia"bles are 
known at all points of network on the row t^ i.e, 
time step t^ and that it is desired to advance the 
computation to time step 

t^+l = + ,iS.t 

Here, the partial derivative at central point. 


M of the network is defined as follows; 
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a(M) = 1 (cej + + af ' + ) 


.D+1 


.3+1 


6 X 


(“) = i>i+i + 4Vi )'(“! + “f ' )] 


(M) 
5 t 


1 ~ Re .. 

2 L\ti ^ 


+ a3+:l ) 

1 1+1 ^ 


+ “1+1^1 


where cx may represent Y, y, A, or B. In this 

conversion, the continuity equation now transforms to 


1 I 


+ yf - 

(yRi 

+ yi 

2 Zht 

1- - i+1 





+ yj )i + 

1 

(A/B) 




2 h X 

0 

-(Tj + 

yh'' )| = 

0 


and momentum 

equat ion 

is, 


g 

2 il\x 

'N.1 


(yj + 

y^^ 
^ 1 


i+1/2 


i+1 


i+1 


2 txt 


(vj:| + vf ^ ) 


(vj + + ■' 


2(ihx) 


Y 


3+''/2 1 rirO+1 


i+l/2 P\+1 


1(^-1. + 


-(vi H- vf' i| + I (s3. + H- sg1 + sgl, ) 


£_ (z? . z? ) = 0 

/^X ^ 1 l+V 
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where 


_ 1 , i -i 

i+ 1/2 - 4 + V^+'' + 7?+1 

1 + *1 ^ 


^''/^h+ 1/2 = ~ jj 73 /A 3 ) ^ ^ (v 3 '+VaJ+'' 


-f- 


(Y?+j /J+1 V! 

^'i+i /‘-ii+i ;] 


=a 


+ 13 + 1 . +1 : 

B- -R 0+1 ^ — 1 + 1 ) 

- 1 B. B^D+1 


( A / B )^’'^'’/’2 ^ 2 

i+1/2 "4 




^+ 1 '"' 1+1 


+( 4 '"Vb 3 -^i j 



J'or frictio 
^sed in the 

superscript 

unknown. 


u slope at 


any point, Manning's formula 


aboi’^e formulae all 
axG unknown while 


the variables with 
those with j+i are 


is 

0 
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These equations constitute a system of two nonlinear 
algebraic equations in four unknovms. By themselves 
they are not sufficient to evaluate the unknowns at 
points i, and i+1,j+1. 

Considering n points on any row, there v/ill be 
n>-1 rectangular grids and n-l interior points. The 
combination of all the rectangular grids provides 
2(n'-1) equations for 2n unknowns. The two more equations 
are provided by the u/s and d/s boundary. 

Now, the channel reach is simulated by i'I*-ordinate 
lines on the x-t planes, a system of 2N nonlinear 
algebraic equations for the solutions of 2N unknowns is 
obtained. They are as follows; 

(y/-i » Y2I2 ^ ^ 0 

•W. IB. 

^'■'.•■'1 ^ n ’ ^ n ) ~ ^ 

'■ n-1 ^^n-1 ’ ''^n-1 ’ ^n "'^n^ 

Y (Vrrc W ) = 0 


0 
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G-eneral Eq_uatlons 

In the system of equations, the values of the 
variables at time step t^ are knom and may be treated 
as constants. The unknowns are denoted by superscript 

(o+D* 

The upstream boundary condition is 

y-i = A y-i -A = 0 

ki “ ° 

The downstream boundary condition is 

° 

For intermediate stations, equation of continuity 
takes the form. 


Si [yo- yi+i> h+]] = (n+i + a) + a 


+ 


i 

4 Ax 


(n+i " yf ('i + h+i 


+d + 


1 ^ 
4 






T. + y. + e)J + ^ 


1 Is t 


A X 


h + m V^- + p 




0 
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and the equation of motion is 


&ijyi, Vi, Vj.+7 


=(^1+1 - yi + a') - ~ (y^ + Ti^.1 + b') 


1/2 1 i 2 « V 

+ ° ’a+1 + a Ti - Ti + e ) 


+ ^ (3fj_ + + h') - 2. C.X - z.^^) 

= 0 

Here all non constant terms indicate the time 
step 0 + 1 . 

where a, h, c, d, e, h, p, w, a’, b’, o’, d’, e', h' ar 
given interms of the variables at time step t^. 


■( yi +1 


+ yI ) 


a’ = y, 


H-1 ^ ^i 


Y? + Y? . 
1 1+1 


b' = 


( y^i + ) 


d 

^i+l 


- y? 

1 


c' = 2Y, 


i+1 
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d = ('4i(yLr yi> + yj) d' = -2T,? 


e = T? ^ - Y? 
^ 1+1 1 




h = aI /3| + i-U/Bl, 


h’ = si- + 


m = - ( 


1 . 

BO 

1 


A. 


+ 


'i+1 


B 


i+1 


iJow the equation comprise. of 2n nonlinear algebraic 
equation in 2n unknovms. They are 


y-l - = 0 


B>|(5''-]j 7^» y 2 ^ '^ 2 ) ~ 

0 

^ 2 ^^ 2 ^ ^2.’ ^3’ '^3^ = 

0 


II 

0 




G--1 (y^i* 72’ ^^ 2 ^ ^ ^ 

G'2 (72’ ^2’^3 ’^3^ ~ ^ 
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In fact it is a very useful information that although 
the equations involve 2n unknowns, yet each equation 
contains a maximum of four unknowns. 


Generalised Kewton*s Iteration Method 

The application of generalised Nev^ton's iteration 
method to the system of equations is made by assigning 
trial values to the unknowns. 


it)i ( 


O' ‘-I’ <j~2’ '-‘3’ • • • 




2 V o-i» 2’ '^3 


u '2 > 




o 


2F 


) = 0 



^2N ^ 



2W 



Solutions are obtained by adjusting the trial 
values until residual vanishes or is reduced to a 
tolerable quantity, let us indicate by suffix k the 
kth iterative value of etc. It is 

desired to approximate the values of the variables 
through the (K+l)th cycle. When the values of the 
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variables are substituted into equation >- - -A, the 
right side becomes the residuals. Let the residual be 
represented by Rl. 


then 



0-1 ( 0 -^ j 
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According to the generalised ITewton’s iteration method, 
the values of the variables are related by the following 
equations 
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where all the partial' deviatives are evaluated at the 
kth cycle, 


d o":! 


dCg 


K+1 1, 
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k+1 k 
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Now the above actuation can be solved for d <j- 
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o'"2’ 


<3-cr 


2N 


and values of 


~V 

k+1 


1 


k+1 

cr- 2 ’ 


k+1 

'211 


can be easily found out. 


This procedure is calculated till the residual is 
negligible. 


Soluti o n of the finite difference Equations 

The same procedure can now be illustrated with the 
flood routing problem. 
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T lie e -uat io ns are : 
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tR — r 


k 


Pi (y^, yj, Tj ) = 


„ / k ,.R R irk \ 

&i (I'i- Vj., yj,- Vi ) 


R' 


k 

2, 1 


/ ^ „k k „Rn _ pk 

^N-1 ^^1-1 ’ ^lT-1’ % ^1,1-1 

n r ^ T ^ 

^N-1 ^^-1’ ^R-1’ ^2, N-1 


% (y^» % ) - 


Tlie residuals and tlie partial derivatives 

are related according to tlie generalised iteration 

method by 
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All the paxtial derivatives are evaluated at. the 
hth iteration cycle. The solution of the system of 

V4. -t lr_|. -j 

equation will provide values of y^ , Yf . The 
procedure vjill he repeated many times until the desired 
accuracy is achieved. The values of the variables found 
in the terminal iteration cycle will be taken as the 
values of the variables for the time step (j+O and 
the computations will be advanced to next time step (;i+2). 

The coefficients of the matrix are obtained by 
diffentiating equa.tions with respect to the 
independent variables and are given below. 
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This numerical procedure was applied for routing the 
flood tlirough the channel. 


i+1 
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^ 5 o f Programme 

In the hook Hydraulic Transients by Streeter and 
Wyle, in the Ghapter open channel flov'/ a problem of 
flood routing has been solved. The same problem was 
taken for checking the correctness of the program. The 
results obtained by the implicit method were compared 
with the solution obtained in the book. It was found 
that the depth compared very well in the two solutions. 
Velocity in the implicit scheme was higher than the 
characteristic method. But the difference was limited 
to the second decimal places. This validates the 
computer program. 



CHAPTER III 


LITERATURE REYIE?f 


The subsidenoe of v/ave in the channel with initial 
steady condition is the phenomenon of prime interest for 
any practicising__^ engineer . The subsidence may be 
controlled by the resistarice and acceleration terms in 
the dynamic equations of motion o?o by a much more 
simple mechanism of pondage in lakes through which the 
flood passes. But before going to the fundamentals 
and detailed anal .ys is, it would be useful to have an 
idea of hovi/ other analysers have explained this 
phenomenon. 

The simplest wave form as given by Seddon l_2lj or 
independently by Kleitz, is the monoclinal rising 
wave in a uniform channel. Scuh a wave consists of an 
initial steady flow, a period of uniformly increasing 
flow, and a continuing steady flow at the higher rate. 
Superimposing on this wave system a velocity v equal 
and opposite to the wave celerity, c causes the wave to 
become stationary and a steady flow known as overrun 
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is = (c - v^)A 


or Q >- cl = constant 


or 


aq _ 

ri 


c = 


i ^ 

B dy 


This is called Seddon law after the man vAio demonstrated 
its validity on the Mississippi River [__2lJ . This 
Seddon’ s law is based on the assumption that q is a 
function of y alone. Secondly, it assumes that the 
wave form doss not subside or diminish in shape which 
is contrary to the experience. To overcome these defects 
concept of kinematic wave v/as introduced by high thill 
and Yfriitham [19]. These investigators defined a 
kinematic wave as one in which Q is a function of y alone 
this imples and that the other slope terms 

are negligible. It can bo easily shown that there can 
be a true wave motion ?men Q is a function of y alone, 
by considering the equation continuity i.e. 
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dr 
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B dy 
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T liGrefore , for an observer moving v;ith a velocity c, y 
vjlll appear to bo constant and so Q will also appear to 
be constanit. This leads to a very important conclusion 
that the Kinematic wa.ves move with the same velocity as 
given by the monoclinal wave. At this stage dynamic 
wave can be introduced. This wave is obtained from 
the dynamic momentum aquation and continuity equation. 

The kinematic waves lead to only one waw?. velocity 

far 

while/the dynamic waves there are two possible wave 
speeds physically this maans that a dynamic disturbance 
i.e. one in which all slope terms arc important v^ill 
propagate in both upstrea,m and downstream directions. 

A kinematic disturbance wrill propagate in the downstream 
direction. 

Actually there is possibility that both the kinds 
of w;avG motion are present in anjr natural flood wave. 

S ince the bed slope term dominates the other slope 
terms, even if the latter ones are negligible, the 
main bulk of the flood moves substantially as a . 
kinematic wave, although its character is modified by 
the other slope terms. Despite these arguments, the 
speed of the main flood wave may be expected to approximate 
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to that of the kinematic %vave and the same has heen 
concluded by Lighthill and Whitha.m. 

Henderson il4j in his papers vjhite discussing the 
flood wave characteristics has shov/n that the leading 
dyna.raic wave acts as a fore runner of the main wave, 
and if P'O > 2, this fore'-runner will bring about appreciable 
rise of the water level before the arrival of the main 
wave. But in most of the na,tural flood wavs problems, 

B o is seldom greater than 2 and so fore runu^r does not 
play any singnificant role. In fact the discussion of 
the fore-runner is of limited practical limited utility 
while de-aling with the natural flow; problems. Henderson 
licis shovm by qualitative analysis that there will be 
appreciable subsidence of the flood wave. The space 
rate of subsidence is itself a small fraction of the 
bed slope. In the most of the flood routing problems, 
to avoid time consuming calculations, a much simplified 
approach w;as adopted by the practicising engineers. The 
flood w;a.ve movement can be simulated by an electric 
analogy in w/hlch a pulse is fed into the system. This 
pulse can be attenuated either by a capacitance or by a 



resistance. When effect of capacitance is dominant or 
hydraulically, when effect of storage is important it is 
called storage routing. Thomas l_28j used graphical 
method of integration of the governing equations namely 
continuity and momontuni equation in his analysis. He 
used the simplified form of equation of .motion by 
neglecting second order terms. This method requires the 
contour map of tlio area and is also crippled with "tedious 
time consuming graphical approach. 

The best known approxinis,to methods for the 
integration of the equations of unsteady flow are the 
so called storage routing 1 20_] methods. In a typical 
raethod the equation of continuity without lateral flow is 
represented in finite difference form as 


£VQ ^ ^ 
1 \ ,1 1 


or 



A, A , is nothing but the change in the volume 

or the change in storage i.e, 

AS = AA . ax 


where S = storage 
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This eqtiation is the equation of continuity. Another 
widelj^ usGd storage routing raothod is the Biiskingum 
method j,20j , This is called Muskingum method hecause 
it Vv’cis applied in the river Iviiskingum at Ohio. 

In this method the storage-discharge relationship is given 
by an empirica,! equation the coefficients of r/nich a.re 
obtained from historical records. This method is described 
in any of the text books 20. ve' ^ . Brakensiek and 
Oomer (1966) gave an approximate approach in which the 
aquation of continuity is retained a,nd the momentum 
equation is simplified to (5y/5x) =3^-3^ by neglect- 
ing the accel8r:.i.ting terms. 

Approxime^te m- thods whether numerical or 
analytical provide solutions which are useful and 
a,ccurate the serious disadvantage is as regards 

the uncertainty in accuracy of solution. YJith the use of 
high speed computers, attempts were made to solve the 
complete equations numorically. The one such attempt 
wfo.s made by M. Amoin in 1966 [3j , Freed ^1973^ Baltzer 
etc. Theoa m^^thods were computer based and so the 
degree of accuracy co.n be taken to any dosired degree 'of 
accura.cy. Amein I 4 ] in report has shown how successfully 
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tho implicit scheme can be used for an irregular channel. 
Out of all the methods namely explicit, characteristics 
and implicits the last one has betui found to be most 
stable. 

The pa,rairietrio study of flood wave propagation in 
a part full pipe wsis done by Peter Ackers and A.J.M. 
Harrison jl] in 1964«. Ihe work was published in the 
Proceedings of the Institution of Oivil Engineers 1964? 
Yol. 28. ThGii’ study was for a pipe flov\7, studiee^ 

only the subsidence of peak depth. In 1969, Mozayaney 
and Song also presented a parametric study on thvC flood 
propagation. 'Thej?- used a sinusoidal hydrograph. The 
wave aRiplitude and duration etc. wore very small so that 
they .are most suited for a small flume. 


In a very recent paper, Sridhara,ii and Mohan Kumar 
studied the effect of various parameters in the case of 
open channels . THe^ha^^enondimcnsionalised the equations of 
continuity and equation of momentum, "flflis study is a, 
good o.pproach towfards the parametcric studies. But 
-bhe^cov rsia and thus 


complicatiiT-g the study 



OH&PTSR. rv 


PORMIATION OP THE PROBLEM 


4-. 1 Oeiieral 

It is obvious from the theory that it is very 
difficult to solve the governing equations analytically 
for various practical problems. The numerical techniques 
are very complicated and for accuracy high speed computer 
is required. Keeping in view all these things it was 
considered desirable to carry out a systematic dimensional 
analysis and determine the variation of the important 
dim,ensionless parameters relating to flood propagation as 
the pertinent independent parameters are varied systematically. 

4.2 Dimensional Analysis 
4.2.1 (reneral 

One of the important dependent variables of 
interest in flood routing through a river is the subsidence 
or decay of the wave amplitude with distance. This subsidence 
is defined as the depth at any section minus initial steady 
depth, non-dimensionalised by the net increase in 
stage at the initial station. In other words this term 
gives the wave amplitude at any section normalised by the 
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initial wave amplitude at x = 0. Mathematically, 


= Relative wave amplitude 

(Ypt = 1 - '^0 




The other important variable is the time of tra.vel of the 
peak, The time lapse between the occurrence of peak 
stage at the initial station and at any stations (downstream) 
is equal to the time of travel of the peak. The time of 
pea]f depth, expressed as Tp is non-dimensionalised by 
x/(V^+Cq) which is the theoretical time taken by the 
wave to reach at particular location. This non-dimensionalised 
time presents an idea of how the wave is moving in time 
space plane. Too large value of it means that small 
amplitude waves are coming first and peak of the hydrograph 
has not been reached. Similarly if it becomes small, it 
mearsthe wave is moving with the peak depth. 


4,2,2 Dimensional Analysis 


The independent variables are the initial uniform 
flow depth Yq, distance of the section under consideration x, 
Manning's n, initial Droude No Fq, time of peak of the 
hydrograph tp, time to centre of gravity of the hydrograph 
t , peak depth of the hydrograph y^, rate of rise of 
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the rising limb described as y^/t , nature of the 

^ ir 

hydro graph. 


Mat he mat i c al ly , 


Y* = f-, (x,yQ,n,y ,Fo,t ,t , 


P’ g’ t. 


shape of the 


hydrograph) 


"P _ ^ / 


^/(Vo+Go) 


hydro graph) 


Dimensionless grouping of the above variables leads to 
the following two equations 


= (— , n, 


'p yp ^“^o 


° yo tp(v^ + Cq) 

shape of the hydrograph) 


(Vo + Go ) 


r (— n 1? P iE_ 

tpW ’ "^0’ + -sr ’ 

^ yo ° ^O 

shape of the hydrograph) 


yp-^o 


where Pq = Froude no. at initial steady flow 


; g is the acceleration due to gravity 

2 

= 9.81 m/sec. 


q„ is initial steady discharge 


y -Jq _ 

-B dimensionless rate of rise of the hydrograph. 
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4- . 3 Parametric Study 

The important factors which control the movement 
of flood waves in the regular channel can be classified 
into 3 groups. 

(i) Channel Geometry 

(ii) Plow Condition 

(iii) Initial Conditions 

The variables in each are as follows 

1 . Channel Geometry 

(a) Cross-sectional shape 

(b) Channel bed slope 

(c ) Roughness of the channel expressed as 
IVIanning's coefficient n 

2, Plow Condition 

(a) Upstream boundary condition- given by the 
inflow hydrograph 

(b) Base flpw 

(c) Peak depth of flow 

(d) Time to reach peak of the hydrograph 

(e) Time to the centre of gravity of the hydrograph. 

3. Initial Condition 

(a) Steady uniform flow 

(b) Steady nonuniform flow. 
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4. 4 Range of Variables 
4.4.'! Eroude I^amber F q 

In most of the natural rivers the flow is in 
subcritical state and leaving the initial reach at its 
mouth , the river flows in tranquil state. 

In their recent paper of A3CB jSSl , the authors 
have concluded that the Froude number does not have 
significant effect on the subsidence in the case of wide 
channels. Keeping in view the large number of computer 
experiments it was decided to use the most suitable 
single value for Froude Number, This value has been taRen 
as equal to 0.3. 

4.4.2 Cross-section_ area 

The simplest one is rectangular shape. So 
rectangular section was adopted. 

4.4.3 Manning ' s n 

Manning's n was varied from 0.02 to 0.03. 14: was 

not possible to go beyond n =0.03 for J^/Jq ^ 3 as there 
was local supercritical flow, 

4.4.4 W ave ampl i t ud e s 

In, most of the flood problems the maximum 
depth of flow comes out to be nearly 3 times the steady 
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depth. So Jp/y^ was varied from 1.5 ho 3. 



The Voorst flood has a value of 1/3 and the other 
one was 2/3. 

4.4.6 Inflow hydro graph 

Inflow hydrograph used are of three types: 

(1 ) Sine curve 

( 2 ) i-i-jely skew hydrograph of type III 

(3) Triangular hydrograph adjusted in ( 2 ) above. 

Sine curve was not used in the further analysis as it 
does not resemble the real flood problem. The suitability of 
-hvely skew hydrograph bound in the direction of zero has 
been discussed with details in Chapter V. 


The eq_uation of the hydrograph used was 


y 


70 


= 1 + 


7 . 
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'P 


where all terms have been defined earlier. 



CHAPTER y 


ANALYSIS OP RESULTS 

5 • 1 G-eneral 

Based upon the dimensional analysis all the 
variables were non-dimensional is ed. Yarious graphs were 

plotted betv/een relative wave amplitude and distance with 
one more variable keeping all others constant. This 
sort of analysis will give an insight into the compli- 
cated phenomenon of wave propagation. 

But the first question arises as to what the 
shape of the inflow hydrograph should be. Mozayaney 
and Song in their paper have used sinusoidal 

hydrograph. But in actual field observations no 
hydrograpli is strictly sine curve, rather they have 
steeply rising limb and comparatively flatly falling 
limb. So it was considered more justified to use 
skewed (to the left) hydrograph in place of symmetrical 
ones. That the hydro graph should be +vely skewed 
one can be shown from the argument which follows here. 
Refer to Pig. 5.1. The hydrograph at the of the 

river is nearly symmetrical. But the falling limb 
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has a velocity = T -i- G.jvhils the rising liml has ¥-3 
where 0 is '’celerity = which of-course will 

die out as it travels upstream. But if temporal 
distrihution of flood wave is considered, it can "be 
concluded that the rising limh will remain steep 
hut falling limb will loose its steepness very quickly. 
Modification of wave with distance and time is shown in 
figure 5.2 which also supports the statements made 
here , 


Before proceeding to further a,nalysis some terms 
are defined here. 


(a) 


R.ate of subsidence is described as 


Relative v;ave amplitude Y. = 


^^peak^x=l ^0 


(y. 


peak) 


x ==0 


^ y. 


The relative wave anplitude gives the 

maximum wave amplitude at any section normalised 

by the same at x = 0, 


(h) Rela.tive discharge rise 
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Q(q ) 

^^^max^xpi 




■max^x=0 


% 

% 
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(c) 


T = 


VCT„ + 0„) 


i.e. time taken by wave pe?,k to reach any 
station since the starting of the hydrograph 
non>-dimensionalised by the theoretical time taken 
to reach that section. 


^ ^ Analys i s of x-ate of subsidence; - 

To analyse this following graphs v^ore plotted 

(1) X with different rate of rise of the 
h 3 ^ dr o graph. 

(2) Y* vs X v^ith different Y^. keeping tp/tg = constant. 

(3) Y* vs X with different tp/tg v.alues keeping 
Y,^ = constant 

(4) Y* Ys X with n. 

5.2.1 Effect of n: 

Plots of relative wave amplitude Y^ and X v/ith 
{Planning’s n as the third parameter indicate that 
tiannings n has very significant effect on the subsidence. 
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Erom the graphs, it ca.n be seen that an increase in 

■/mannings roughness coefficient result in the increase 

w 

of relative wave amplitude. 


Eor n = 0.02, = 0.550 for tp/tg = 2/3 at x = 16 

for n = 0.05, Y* = 0.685 for tp/tg = 2/3 at x = 16. 


Therefore, the effect of Mannings n is to de.creass the 
subsidence. When Mannings n increases the wave pa,ss 
without getting subsided. This may be 
explained from the momentum equation. Energy slope 
is directly dependent upon mannings n and so with the 
increase in manning’s n, energy slope also increases. 


5y 

ST 


+ 


Y + 

5X 


o- 


g Sq - g 


how, the right ha.nd term dimini-sh and so do the left 
hand terms, i.e., the acceleration term decrea-ses. This 
appears to be the possible explanation of this effect 
of ma-nnings n. 


5.2.2 Effect of Wave Amplitude 

The variation of wave amplitude Y* with X 
keeping tp/tg constant is plotted in the figure 5.3 a.nd 5.5. 
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Erom these ciirvGS it is clear tli?at relative v/ave 
amplitude increases with increase in Y„, . This 

vv 

increase is small a,s compared to the effect of 
Manningb n. 

Secondly, it can be seen that the steepness of 
the plotted curve decreases with the increase in distance 
X. The wave in later reaches tend to become flatter, 
and the subsidence in the later roaches has been found 
to be more. In the initial reaches Y* and X vary in a 
straight line pattern. The same conclusion has been 
confirmed by Sridharan [26 1 in his p-aper. In fact this 
straight liness is true for small length of the river, 
Lyer used 30 km of channel length and his conclusions 
also support the straight line nature, 

5.2,3 Effect of tp/tg: 

The plots of Y^ with X for different tp/tg 
values are presented Eig. 5.6 and 5.7,Y* decreases 
linearly with X in the initial reaches. The line tends 
to become flatter vrlth tncrea,se in dista,nce. The rate 
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of chringe in slope decreases with distance. The 

effect of tp/tg is to increase the relative wave 

amplitude Y*. This effect is more pronounced than the 

effect of wave amplitude Y,„. 

w 

Suhsidence is more in the later roa-ches. This 
can he e3q5lained as foLlo?ra. When the wave has 
flattened sufficiently in the lower reaches, the waves 
become short duration ones and rate of subsidence 
increases significantly. 


5.2,4- Effect of Rate of Rise of the Hydrograph; 

Rate of rise of the hydrograph is defined as yw/tp. 
■from the computed results of rate of rise of hydrograph 
it V'jas found that this term practically has no effect 
on the relative waA^'e amplitude, fig. 5. -iiiclearly shows 
this. In fact this rate of rise is taken into account 
jointly by 3 rw/yo and tp/tg and so this term can be 
ignored. 

5.2.5 Effect of shape cf the Hydro graph; 

Keeping tp/tg, yw/yo, n same various runs of 
the computer experiment were made by using triangular 
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hydro gr "ip lis . I'rom the computed results it was found 

that the shape does not play much role in subsidence 
phenomena, il though in case of tp/tg = l/3, the 
triangular shape showed a liutle more subsidence than 
the +vely skewed hydrograph,' ' This can be explcuined as 
follows In the case of triangular hydrographs of very 
steep rate of rise, the falling part gives more depth 
than the +vg skewed hydrographs and so the cumulative 
effect of the triangular hydrograph is slightly larger 
subsidence, Thisyshovm in figure 5. .'O- The above 
mentioned two conclusions are very significant. Now 
only n, yw/yo, tp/tg remain to be studied. 


5.3 Rating Curves; 

A graphical representation of stage vs discharge 
is called a .rating curve. Such a rating curve for a 
typical ca.se is plotted in Pig. 5,3iCfyThe rating 
curve clearly indicate that there are two values .of 
stage for one discharge. On the rising side the 
discharge is more for the same stage as compared to 
the falling side. Secondly with increase in distance 
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looping curve tends to "become thinner meaning there'by 
unique value of thie discharge for one value of the 
stage. (The first observation leads to the belief that 
tho discharge is not a single ■val.ued function of stage 
or discharge is a function of some thing else other 
than tho stage. In fact during the period of rise of 
the hydrograph discharge is more but after the peak 
of tho hydrograph, discharge falls down owing to the 
reverse slope of water surface. 

5.4- Analysis of Time of Occurrence of Peak Depth; 

lime of occurrence of peak is defined as the time 
of occurrence of flood peak at a particular section in 
the reach minus the time of occurrence of flood peak at the 
initial station. This is the actual time of travel of the 
peak, 

5,4.1 Effect of n: 

Prom the computation of results it became clear 
that Manning’s n plays practically no part in the time 
of occurrence of peak. Pceviously in the analysts of 
relative wave magnitude this very term was found to be 
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a dominating factor. It confirms that in the left 
hand terms the effect of temporal 

acceleration is little as compared to the spatial 
acceleration. 

lime to peak is non dimensional is ed by the 
theoretical time taken by the kinematic wave to reach the 

section and the plot indicates that xy^Yo 4 - ' ^ ' ) 

independent of mannings n. Therefore it can be 
concluded. that time to peak bears a constant ratio with 
the theoretical time taken by the wave to reach at the 
same section, 

5.4.2 Effect of 

■ye Y 

Graphical representation of 
with third variable indicates the following 

points. One is that in the initial reaches there is very 
steep fall in time to peak. This is probably due . 
to the discontinuity of the term Tp/k|[Yo+Go)) at x = 0. 
Secondly, in the later reaches the ratio begin to 
decrease with distance or time to peak is quickly 
reached at the section under consideration. Thirdly, 
the Tp/^Yo + C5o|) value tend to become constant in the 
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further reaches, fourthly, the effect of Initial 
wave ratio is that with tho increase in Y^, 


Tp 

3vrynTo) 


increases or lag of time to peak v</ith theoretical 
time increases. So it can he concluded at this stage 
that the effect of initial wave ratio is very significant 
in the analysis of tinn to peak. This property can he 
of great use in the prediction of time to peak as will 
he shovm in tho next article. 


5.4.3 Effect of tp/tg; 

If time to peak Tp is nondimensionalised hy 
theoretical time r8q.uired to travel a distance is 
plotted with distance x, it is found that tp/tg of the 
hydragraph is not very important variable. 

Therefore it appecors that the wave travels 
v/ith a velocity which is independent of the tp/tg 
ratio of the hydrograph. This conclusion has been 
confirmed hy Sridharan [2o3 also. 
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5.5 Prediction of Ip and Tp 


5 .5 .“I Q-eneral 

In any flood routing problem, the engineer is 
interested in finding the magnitude of the meocimum stage, 
maximum flood and time of their occurrence. By the computer 
solution based on numerical methods, one gets detailed 
picture of the surface profile, stage at any time at any 
section etc. which are not of much use to a practising 
engineer, further the computer programmes are very 
complicated and require a lot input data in an specified 
form. All these factors call for a direct short-cut 
method for the prediction of peak stage and time of its 
occurrence . 

After the detailed study of wave propagation based 
on dimensional analysis it is now possible to drop one or 
even more of the variables as their relative effect on the 
dependent ^’•ariables is less. 

5.5.2 Principle 

from the previous study in this chapter, it can be 
concluded that for the rate of subsidence of relative 
wave amplitude is not dependent upon shape of the h 3 rdrograph. 
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rate of rise of the hydrograph, wave amplitude Y^. fhe 
most important variables affecting Y* are Manning's n 
and tp/tg* A. plo 0 of i "“Yit* on natural scale and Y on 
logarithmic scale indicates that after some initial 
points, the points lie approximately on a straight line. 
This plot indicates that 

1-Y* = Y+ = c + m log^Q(X) 

where c and m are functions of n and ti^/t^. 

P & 

t 

i.e, c,m =f-i p(n, r-^ — ) 

tg 

This is a very important conclusion and forms the basis 
of prediction of Y* , Similarly from the analysis of 
time to peak stage, it is clear that it depends only 
upon Y^ and X. This forms the basis of prediction of 
time of peak stage. 

5.5.3 Procedure for graphical prediction 
Rgla-tive wave amplitud e 

On a semilog paper, X is plotted on logarithmic 
scale as abscissae and 1 -Y* = Y + as ordinate on 


arithmatic scale. 
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After some of the initial points, which correspond 
to the straight lime portion of the Y* vs n plot, a 
straight line is approximately fitted through them. 

The slope and intercept of the lines can be calculated 
by taking any two points on the straight line and 
solving them for m and c, where m = slope of the line, 
c = intercept on axis. 

ID \ 

now, m, c = f-| 2 ^ ^ 

* t 

g 

So, m and c are plotted on a simple graph paper with n 
on abscissae and tp/tg as the third variable. 


Since soma error is caused due tc neglecting 
of effect of "fepAg ^ suitable additive constant should 
be used for each y^ value, 

T 

(b) Time of peak T^<- = — ^ 




A plot of ^pAg 0^1 ordinate and X on abscissae with 
y^ as the third parameter is prepared. The dotted lines 
on the plot indicate the interpolated variation Tp/tg and 
X with other y^ values. 

Additive Constant for Y-)f 
V, -r- 0.5 1 1.5 

i VV 

Constant to __ 
be added to ^ 


0.015 0.02 


0.035 
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5.5.4 Steps 
Data required: 

Initial steady discliarge , Manning's n, t„/t 
of the inflow hydrograpb.. Slope of the channel s 
distance of the section under consideration x time to 
centre of gravity of the hydrograph, 

(1) From Manning's equation initial steady depth 
should he computed: 

I 1o ^ 1 5/5 

i.e. y = 1 1 

S 

X 

(2) For the given x, X = is calculated. 

yo 

• £. - - 2 . 

( 3 ) Using the plot of 2'ig. ■ , , , against tp/ig and 

n value m and c are calculated. 

( 4 ) How, 

Y+ = m log-jQ X - C 

y+ can be calculated. 

(5) = 1 - h 

( 6 ) Y* = Y* + additive constant. 

This additive constant has to be interpolated from 
the table given above. 
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^ypeak^x=i ~ 

(7) Y* = 

^^peak^x=0 " 

from this (ypeak'^x=i calculated. 

(8) Now coming to the Fig. S3. 

From this plot for the given X and of the 
hydrograph, ^p/tg can be noted down and tg being 
known can be calculated. 

Illu stration 

For X = 500 Kj^ , n = 0.02 , Y^ = 1.0 
500x1 0^ ^ 

^ _ = 10^ t- = 12 hrs. 

5 S 

. ^p 2 
for n =0.02 and 

h 3 

m = 0 .044 
c = 1 .845 

1+ = -1.845 + 0.44 log-io ''0^ = 0*355 

= 1-0.355 = 0.645 

Y* = 0.44 + 0.02 = 0.665 

y = 0.665 X 5 + 5 = 8.325 

which resembles very closely with the value of 8,26 m. 
For time to peak, Tp = 12 x 3.30 = 39.6 hrs. 40 hrs . 

From the computed results also, we get Tp = 40 hrs. 
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Limitation 

(1 ) The magnitude of numerical damping is not known. 

(2) This is strictly for skewed hydrographs but can 

be used safely with triangular hydrographs 
for "tpAg upto 1/3. 

There should not be local supercritical flow. 


(5) 



CHAPTER VI 


6.1 Conclusions 

The parametric study on the propagation of a 
flood ware in open channel leads to the following 
conclusions:- 

1 . There is significant subsidence in open 
channels. 

2. Bflanning's roughness coefficient has 
significant effect on the subsidence expressed 
as relative wave amplitude. 

3. The non-dimensional time to peak is found 
to be unaffected by the manning's n. 

4-. Shape of the hydrograph does not matter much 
in wave subsidence. 

5. Rate of rise of the hydrograph does not 
affect the rate of subsidence. 

6. Inflow wave amplitude has some influence 
on subsidence, 

7. Wave peak duration tp/tg has pronounced 
effect on subsidence. 

, Hon-dimensional travel time is very much 
affected by the wave amplitude and is 
not affected by tp/tg. 


8 
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FIG-5-2 SPATIAL PICTURE OF SUBSIDENCE WITH TIME 
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APPENDIX 


PLOY/ CHARTS 




gLQW CHART 


Main Program 


(start 


^ead 

M, lT,L,dx,dT,TPIN,SO,&,Tzero 
Xzer 0 , tg , tp , tp , Y q , dishi 


z 


stage- (I)=f(t) 
1=1, MN 


Jl 


. , . dishl 

Ylir(l) 

S ^(yo) 
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Ey(l)=YIN(l) 
I =I,N 

E 


T(K+1 )=G}(K)+dT 



I Call 

Inter • 

I 

1 -.. 

} C^l Dwns^ 


T 

jL 


j Set Coefficients of the! 
\ matrix and Residual ^ 

Vectors ; 


Jl 


Call TCTVBR 


I 


Compute Error in Trial 
Values and Computer Values 

^ j ^ - 



Type depth, \ 

time, discharge) 
velocity ! 


\ 

L 

Calculate new i 

trial values of 1 

Vel. and depth 1 

for next time steps 

> 

/ 

' Redefine vel, 
de-Dth . 

and 



\ 

\ Time 

= time+ 3>t 

























Subroutine IlOTER 


Subroutine lUVIR 


Calculation of vel.and 
depth at all distance 
points 


Redefine depth and vel, 
at all 4 points of a 
grid system for all the 
distance steps 


DO-I = 2 , 


Call’- 

/a/ 


Calculate i.,B and n at 
all distance steps 


Set the coeff. 
of the matrix 


Calculate the average 
values of n and the 
first and 2nd derivati- 
ves of A ahd B 


Subroutine ’Di-- 


Calculate the coeff s,. 
of the partial deri- 
vatives 


Calculation of vel. and 
depth at downstream 
point 


Calculate the values of 
the coeff. at the Kth 
cycle 


Compute the coefficients 
of the partial derivatives 


Calculate the values of 
the residual vectors 


Set the coeff, 
of the matrix. 


Repeat computation 
for all 1,1=2,111 


Return 
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PVCiO =X(i'L) 

iiH/vyKixRhyti 

?F(irKAy-50)330,330,335 

T y p t 3 3 J , n H A 1 . t V 1 1 5 , e y t f 1 , fc y 1 1 ) , fe. V ( i : j 
FJPMAitlh ,i?»Ab-;ll .3) 

CALI. tXiT 
CilNllNtlE. 

03 240 I 

FKROR =(EVClD-'Vi(l))/fiV(T)»l()0,0 

3f (MBS (ERROF) -2.01 24o,'24J,b5 
CJKIliUiiS 

Ip(lP-TARii) 245,245/255 

CDMIwUe 

TP=IPt2.0 

lypfc 3i,n?4. 

TiPfc; 23 

FaFvAl(5X,F7.2) 

03 25w 1 =l,L 
fcA(i)=AH(tyu))' 

Fy{I)=EAlil»FV(l) 

TiPc. 25,t)lSi'(i),£i(i),£\?a) 

03 260 1 sl,0 

C3Nn-M!£ 

rXRAYs:! 

DYcii sEyaj-yi-ia) 

Ov'tiJ sEVCll-VlNfl) 

OJ 27U I =l,w 
Y/Ul) sEUI) 
yiNii) savd) 

TlBe =iIBE4DfcM 

IFC2-'lIi4£) 300,2^0,280 
03 285 i *1,0 , 

Eidj =c.yc.r)4iiyt jjt.5 

fiVCl) . =EV(IJ+0VCI)»,5 

F3»^AiC35A, 'SfRe A'-iTLOx KOUTiMU Bi InPoICIl i4ElHODV33X, 
211 1 If/ J1 , . 

F3HLAX(bX,ll(2X,F9,4)) ^ 

forma'icisx/valuss of the variables*! _ 

FaHUA’ltbX, *LrSTAKCE' ,3X, 'DEPTH*, 5X, 'DISCHARGE*) 
FJHLAlCbK,F6.l ,4K,F5,2,4X,F10.3) 

CALL EXIT 




SoBRQUrii'iE UPSTR 


’’cD“’0'v/A9LQK/'X(150)criSXtl50),YINU50J,l(IHU50),Ey(l50).Ey(15 
2JGtl50>,OlSCH{i50),TUb0l,vr(I50),y‘I’C150l,y(I50),Eg(l5O), ! 

^ V ^ * iTv k x' ^ I ■ H MM A ^ « e? % 1?* N ^ ^ C /l 



65 

DO 7 0 I 

si , i«L 


xc U 

=0,0; 


D3 70 3 

= 1,^''L 

70 


=0,0 



71 


7 4 
80 




GCI TO ao 
Ij “"Ij + i 

vSURFL. =STftGe:(u-l) + ( TIME-T(rj-l) )*(iTAGt!;(L)-SIRGf'UL-l J)/On?lE 

SlfCU =SUKFL 

I =1 

CaGL iNVKfi 
AU,U =DftlCli 

AC1,2) =DRIDT 

A (I, 3) SOKIDV 

A(Z,1) =0R20U 

AU,2) =DR2DT 

A(2,3J =DR20V 


T)^P£ 14,0Rll)U,DRlDT,DklOV.DK2OU,DR2QT,DR20V ^ , . 

f JRtaT(2X, 'A(iel) = *Ix,F5.^,lX,'A(l ,2) = %U,F5,2»U,'AClf3)s' » 
lPt>.2,'A(2a) = S5.2,'A(2,2) = *F5.2,2X, 'At2,3)s*FS,2/80(lH-)) 
Ctn = 0 R 1 DIJ»Li + i;iR1DW»w4-DR1D1*R-UR1 


ix, 


1% 


JtU =DRinU»U + URiDv*V4-DRlDl*R-uRl 

Ct2J =DB2DU*U+DR2DV»V40R20T*R-OR2 

TiPi;: 24,C(1),C(2) 

FJRiaT(3X,*C(l ) = \Fi,i,2X, 'C(2)s',Fb.l/20ilH-J) 

HtSTIjftiv 

STD 


Su3R0uTlMfc INTER 


150) ,l>lST(i5u) .YlNtlbO j f VlNU50),Eva50),EX(l50 
ajGtIbOlrOlSCHdSo) ,TC150),vrclb0) .3fTCli>0 3,y(lSa),Ea(150)f 


, i k. j r i J i 

ijytlbO),uV(tbO) ,iiTAGt(l5u) ,Krt(l5o) 
CJNMO.»/i3BiiOr;/G,JTX.'4E,SC,l>TPX,DXi)T,DX,l3T,TX‘4efPRl j.yR2#5,R» 


i- J" nui'i / tjrjijurv/ <3 tiJ 1 Aflt, » OV. r i- *iv A f UAU 1 f i/A » w -W * A *>1:. » UA j, » wn* 

? Mi , T R 1 Du .DC 1 D V , DP. ins , uR 1 D 1' > UR2Dli , UR2 D , aR2 OS , DP.2DT 
COD.iaw/C:aLOK/R*i,N,l.,Nb,MN,I,'<M ' 

C.3MMDu/L)RLDK/sa,F,LFvr!,FDF:V'H,XFIu,i)EDX,DEM.XZERIj,0FIR 
r';i,viMU.^y(.''Wi.flK/Af0hh.9ftS3 .f*i f 2653 .C(2b5) .Xf 2653 


J 1 

Igo 


. r*Lir, v»J# r ur.Vii # km #uc*ija^ u&ua # a^ii 

:u.v,Mt)i*/FBDyk/A(2fcb,26b) ,Ci (265) ,C 1265) rX(26b) 
fXPE il 

FjKy.AT(3X, 'C(-'.),C(M+1) ,ACM,J)#AtM,J + l) ,A(»,J+2),A(N,J + 3),,^ 
.U((.i+i,J),Atri 43 ,J+ 1 ),A(;v,+ 1 -J+ 2 ),A(; 3 * + 1 , J +31 ARE AS fOliDOWSV) 

CQ-’PUi ATIOR AT iniERTOR PuIfvTS 


23 


22 

J5t) 


Ou ISO 1 =2»i'iN 

CAJ.L. INVER 
J =2* 1-2 

SDRIC'S 
U^mI+N =uRtOU 
AlN,.J+2) =DR10T 
A('1mJ+ 3) spRlDV 
A(M+1,J) SDR2DS 
A( HI ,.7+1 ) = 0P2PU 
AiD + l m**^2)=DR2DI 

AOHl , J + 3)=DR2UV' . ^ , 

C(?H =i)RlOU»i..+uRlDV*v+DF.iDT*R-i)Rl+DRlOS*S 

Ct HI) =D?.2rDl'U+DP2PV*V+UR2Dl’^K-0R2 + t>R20S*S 
TYPE 23,CtN).C(.Vi+i) 

F l3 I'' T C 2 X 2 F' & 1 / ) 

T 3 :PE 22 ,A? J) ,A(i‘ ,J+l) ,AC«, J+ 2 ) #A(j 1 ,J+ 3 ) ,A(M+l,J) r ACM+l, J+i) , 
lAtfi + 1 , J + 2 ) ,ACM + 1 ,J+ 3 ) 

KPE 11 

FjPMAi(2X,&F5.1) 

ONiiHUE 

RETuRi-^ 

!?;u) 


THIS SUERDUXIoE IS FOR CAeCULATIO^S AT UO^NSTREH bOUNDAKY, 
SuBRtlJTi-'JE Owns 



CjM-iOM/AtiJ^OK/Zt iSlO.niSi tiSu) bO ) , y t (i( 150) , KV ( 150 ) ,K)t C 150 J 

■ HlbO),Y]PC 15 p),YT( 15 n),^(l 50 ),Ea(i 50 ), 


20 G(lb.)),DI.^CHtl 5 .o),VUbO),YTC 150 ),YT( 15 n),^(l 50 ),E<, 
ja y t I 5 P ) , i> V c 1 50 ) , iJT AGb ( 15 0 5 , EA ( 1 50 5 
C'j>*)aK/ 5 HL.UK/(*,oTXMii,;iO,r.'iDX,uxux,nA,uT,'PiMfe;,oPi,i>Pi 



0 P 2 

DP 2 Di/,uR 2 f)S,uRi!r):i’ 


0 E»iX , D t (i i’ , AZE RU , uF t tl 
t?h 5 jrX( 265 ) 


CU IpiirArTUN A1 u/i> tiOuMuAKY 


o 

R 

DADV 
Vi 
Ai 
A 2 

W 

Cu=CltR) 

Dk1=R-5.0 

OKlOTsl.O 

oaiuvso.o 

Al'JLi, yb-UsORlDT 
A(!H,,iyr.j =UR 1 UV 

Ct'U) =DRini*r;+i)RlOV*P-i)Rl 


=EVU.) 

sfcft 

sDAYCR) 

sAPiHj 

=rtl*Rl 

=P*A 1 


RriruRii 

SRD . .. 

: 3.0 S & £ 5 S S,S S $ s £ & S 6 Ssi e $ E i $-6>. S $ £ 3 & $ 5 S &.$ S:.$„S $ S $ $ S S fo 6 $ S $ $'£».S t 


;?. i; -i 


r<i£p ^aUbPUUT liyE CM-CHljAlFS IPE COEFl I C XEN £ MO KESlDPES 

iiljiJrtOUJ’lME l^VEK 


CQiAMOu/ABbOK/Zr 150) .DISK 15U) , VINCI 50), V INC X5U) ,EVC150) ,EY(150 
2 JGtlS')) ,DlSCrifl50),TCJ50),VTCJ50), YTtl50),i!(l50),EQ(150), I 

30YC.150),L>V(15O),&TrtGE(X505 ,eA(l5u5 ^ . 

C jHM’Jw/at3i,OK/0,OTlM£,SO,tVl'DX,L>XL)T,DX,lJT,TIME,DPi ,UP2,b,K, 

2^,0, DRi i)0,DRinV,ORlOS,DRlDr,OR2L>U,i)R2DV,DR2DS,UR2DX 

“ j N: « r: / c e L 6 N/ i, N , f' , L , 'i l , w f , i , 

n. y . •« 1 ni, y tiM iri c'tjn ir’i tr ii'/.? V'c»T.<. 


: a di-i / oBuOK/so , elev n . ri,Fy c; , xfin , oeox , deux , xzero , dein 

- ‘ - (205,265),C1(265),CC265J,X{265) 


Ci;P' 5 Ai>/F 3 tinK/A( 205 ( 


SOLUTIOiv OF SYSTEM OF EWUaTIUNS WIIh 4 NONZERO TERMS ON 
THfc oANDEO Matrix 


iJ 


u :»EVCX) 

V =EV (1 + U 

5 =EYcn, 

R stYCl + n 

TiPi. ♦.U/K 
f iPlu 39 

EJRwAfClOX, 'S= 

3 -i 
Rb 
Uti 
i/js 
Al 
A 2 
Ml 
32 
A,3l 
Ab? 

3il 
332 


K® 


: 3 i 

C32 

Cll 

c.l 


=yipci) 

=YlKCl+i) 

sVIt-Cl) 

=VIMli+X) 

=aRiS) 

=ARCR) 

=BRtS) 

=rtKca)^ 

=aRCSB) 

sARCRB) 

sjiRtSii) 

=BRCRB) . , , 

=CI(3b) + (CN(Si3)-CI(S0))*UGCl) 
=ClCRb) + (CBtRu)-f'I(kB))»DGCl+l) 
=CI(S) 

=C>‘tS) 


Cld 

C-i2 

ZjX 

r.'}2 

h 

'U 

UaIuY 

Dk2D^ 

DiiUir 

0A4l'V 

Vf^ 
n 
02 ■ 

Vl) 

El 

E2 

P.di 

Ki:a 

QJ 

BO 

AO 

SfM 

Sfi 

SF4 

OoFl 

i;'oF2 

3r 

n;i4 

!)di 

0(52 

OBI 

0A4 

(■)A3 

Oh2 

DAt 

OKI 

0B2 

•2 

Drtll'S 

OBlUT 

OKI DU 
0K2UU 
DK2DV 
i)K2DT 
0ii'2DS 

KiiTiJR.'. 

S. ,■{ D 

SUBUOUT 
0 i'-vL?-BT 
EiJR/,AKU 
TOf.sy.O 
FB=0 

. 

1) J o 5 J 
JjfsJ + l 
JJ=wM+A' 
HiaAs.j. 
I fsJJ-J 
l4i=.J + 4 
XFCul-'4 


=C1CPJ 

=C^Hft) 

=cTi + (CRi-cii)*L)Ga) 
sc 12 + tCR 2-C 1 2 ) ♦OG 1 1 + U 
sg(l) 

sDAlf(S) 

=DAy(R) 

=D‘3y l.S) 
si)ex(K) 

sIVDmYiP J 
sDOA/C+J 
sDAi(Sbj 
=OAXf KB) 

5/ + V3-0-UfS 

sv* +u« vB**ntj 

sh-S+KB-SB 

sk+S-5B-RB 

srf+U+VB+OB 

suAlOl 

sl)A2.DY 

sUAilJY 

SDA40Y 

sy 1 ( D /iB 1 +V/ A2+UB/ Aril + Vrt/ hB2 ) 

= i./ei + l,/E 2 + i./Eiil+l,/e;B 2 

sAl/Ei+A2/E2+Aei/eBl+AB2/EB2 

=C01*C01»U4ABS(U)*(S) +*(-4,0/3. 0) 

=Cn2*Cti2*v*AfjS(v)*CH) ♦♦(-4,0/3. 0) 

=C:Bl*CBl»uH*At5{Uni5)*C5B)**(-4.0/3.0) 

=Ca2*CB2*VB+AtoSCVB)+CHB)*+C-4,0/3.0) 

sCUi*COl*U*(S)**(-4,O/3.0) 

=Cn^*CD2*V*K**C-4.0/3,O> 

=BFi+SFA+i>F3+BF4 

=08^0Y/(Ki*ii1 ) 
si.iB4Df/(r2*E2) 

=OB2DY/i52 

=BAiDY/Al 

sDAlOY/iAlfAU 

=OA2DY/(A2*A2) 

s()A2i)y/A2 

= ,25*DTt)X+V0*Dl+O2+DTi)X*.25*Vl*A0-,5*Ql*n<r*BD 

sDrDX*V2/G+,25*^n*Yl/G+,5*r>x*iD/G+01-DX*2*S0+,5*DX*SF 

S.25+DTDX 

= J-W*VP + Vl*CDAlDY/El-Al*Di33)*4 + ,5*Dr*t>B3*Ql 
sw*yC>+l+W+Vi*(r)A2l>y/E2-A2*OB2) + ,5*Or*Ql*D22 
=,2b*DrDX*Dl+uTDX*AD*.25 
=,25*DTDX*Di-uTQX*AD*.25 
5(l)XDr-,5*(U+UH;J+.5*DX*Ql/Al)/G+DX*0SFl 
= (OXD’1 + .5*(V+VB) + ,5 *OX*u 1/A22/G+DA*DSF2 
= l-.5*nx*vl*V*BA2/5+l)X*(4.0/3,0)*SF2*(Dei-DAl)*.5 
= -l-.5*l-X*Oi*U*iJA3/G+uX*(A,U/3,0)*5Ft*(DB4-OA4)*.5 


iMfc, I'lVEH (A,B,n,XS) 
jw A ( 1 ) ,n( i ) 

SOLiUTiOh 


sl,f4 

+ 1 

J 1!>, 10,12 



J u 

'^£U 

3j 

3^ 

35- 

0% ■ 


4^ 


00 30 tsO,H 

SfeAiiCH FOk WAXImUH COEFFICl&Ml' COLiJMi^ 

U = J 1+1 

P' C A;?S ( 0 IG »\ J - ABS ( A ( X J ) ) ) 70 , 30,30 

'3jlGA = .v( 10 J 

1 A=i 

C'i^'TlKijfc 

('c;5i tOrt i>IvriX IiFttS TnA.*. lOGFrtAi^CtCSlAGULAft 'MATRIX 
I F t Aas ( b 1 G A ) - XOi^ ) 3 6 , 35 , 40 

K ijsi ' 

KCiTUR.v 

I+TtRCHANGE RUFS IF ^’liClSSSAKy 
M®0 +(-!’•' (J-2) 

IX=:I?»A»0 
02=0 + 'l 

lFCju2-i'0 42,42,44 
El 2 “■ !< 

03 50 KsJ,l,2 
li = i:i+!. 

T2=I1+IT 
SAVt^Ktll) 

A(Ii)=Atl2) 

/UI2)=SAVE 

LUVlDt; tiQUAlION XiEmDIFG CDEFFICI KN'X 
P 104 

AEXl )=AEIi)/BiGA 
SAVbstCifAX) 

OlIrtAXlsBC J) 

A(J)sbAVE/BiGA 
FlilMluAlE VEXI VAHTEiPLE 
IFEJ-r.) 5b,7l/,5b 
TuS=.\ + (J-l) 

EilsJY + A 

lfCb3-/'.) 57,57,59 

f j 3 s f i 

00 o5 I.<=JY,r.3 
IXJ = J.0S+IX 
T1’=0-1X 

00 60 JX=.TY,G3 

IX, 1XSi;i*CJX-l) + lX 
0JX=IXJX+1T 

A(IXOa)=ACIXJX)-(ACIXJ)»ACJJX)) 

AlIXlsB(lX)-rB(J) + AtTXiJ)) 

3ACK SOL, U 11 UN 

X. f=H + + 2 

U 0 H 0 J = 1 , b Y 
lAsil'-J 

no bu f.si.j 

BE Jb:)=BCIB)-ALlA)»8tIC) 

PO lC=iC-) 

kiSTtiKiv. 

E n 

■ • ■ 

!: Tills FiJivClIO*^ SuBPF0CtliAr.rt«xiTRN^ARBX’Ai^"3iHER VAbTSbUes' 


5 b 
52 


5'i 


59 

3/ 


6w 

fi5 

70 


FJOCTJOiv; AREY) 

Ak=Y 

Rii,TURH 

E.II.D ■ 

FUiVCTiOn BRCYJ 
BK=1 .0 
RiiTUKfii 



Ktincncsfi ci(V) 

ci=o.u2aoo 

RESTURW 
E,4D ■ 

F-’lJMCTiO!'! CNIY) 

C*=U,0200o 

RfiTURN 

g-^U 

rUNCTiON OAY(X) 

OftYsl.O 

RETURN 

Ki^D , , 

function OBYfY) 

OBY=0.0 
RETURN 
. gt<D 

function DDAYCY) 

D0AY=0.0 

RETURN 

il4FLaN HYDROURAPH/iS, ^ . ^ 
iNFliUW^HYDROUPAPH’iS 
""siiilpi=”""”2louoooE-o3‘jo= 


Tl«E= J. 


UISXAJjC 
j • 0 
2. 1-' 
4.U 
o • <J 

t),\J 

1U,U 
1 2.0 
1. i , 0 
I o . it 
1 . 0 
2Q.0 
22.U 

24.0 

2o,q 

A6 9 0 

30.0 


t. 


.00 

VALUES Ot 

depth 

5.7U 
b.bS 
Ij.-iA 
b,3:i 
5.27 
5.21 
*).17 
5.13 
'5.10 
b.OR 
5 , u6 
5.05 
5.04 
5.03 
5.02 
5.00 


lOfc VAHIABL 
uISCm? HGE 

I2,fe^7 
12.17 i 
11,821 
11.544 
11.329 
ll.loO 
U.02? 
10.92?> 
10.847 
10.7S5 
10.733 
10.702 
lO.hTO 
10,670 



